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Abstract 

We study the topological charge distribution of the SU(3) Yang-Mills theory with high 
precision in order to be able to detect deviations from Gaussianity. The computation is 
carried out on the lattice with high statistics Monte Carlo simulations by implementing 
a naive discretization of the topological charge evolved with the Yang-Mills gradient 
flow. This definition is far less demanding than the one suggested from Neuberger’s 
fermions and, as shown in this paper, in the continuum limit its cumulants coincide with 
those of the universal definition appearing in the chiral Ward identities. Thanks to the 
range of lattice volumes and spacings considered, we can extrapolate the results for the 
second and fourth cumulant of the topological charge distribution to the continuum limit 
with confidence by keeping finite volume effects negligible with respect to the statistical 
errors. Our best results for the topological susceptibility is x = 6.67(7) X 10 -4 , where 
to is a standard reference scale, while for the ratio of the forth cumulant over the second 
we obtain R = 0.233(45). The latter is compatible with the expectations from the large 
N c expansion, while it rules out the 0-behavior of the vacuum energy predicted by the 
dilute instanton model. Its large distance from 1 implies that, in the ensemble of gauge 
configurations that dominate the path integral, the fluctuations of the topological charge 
are of quantum non-perturbative nature. 

1 Present address: AEE INTEC, Feldgasse 19, 8020 Gleisdorf, Austria. 



1 Introduction 


The discovery of a fermion operator [T] that satisfies the Ginsparg-Wilson (GW) re¬ 
lation [2] triggered a breakthrough in our understanding of the topological effects in 
Quantum Chromodynamics (QCD) and in the Yang-Mills theory [3J 0J |5] [6, 7j. This 
progress made it possible to give a precise and unambiguous implementation of the 
Witten-Veneziano formula mmmm- 

In lattice QCD a naive definition of the topological charge density needs to be com¬ 
bined with an unambiguous renormalization condition. The cumulants of the charge, for 
instance the susceptibility, require also additional subtractions of short-distance singu¬ 
larities to make them integrable distributions. If the topological charge density is defined 
as suggested by GW fermions HI 131 HI , however, its bare lattice expression and those 
of the corresponding cumulants have finite and unambiguous continuum limits as they 
stand, which in turn satisfy the anomalous chiral Ward identities 13 HUE]- By combining 
a series of those identities, the cumulants can be written as integrated correlation func¬ 
tions of scalar and pseudoscalalar density chains or combination of them 151 El HU- In this 
form a particular regularization is not required anymore to prove that no renormaliza¬ 
tion factor or subtractions of short-distance singularities are required. These expressions 
provide a universal definition of the susceptibility and of the higher cumulants which 
satisfy the anomalous chiral Ward Identities |7J. 

Recently a new definition of the topological charge was found | T2] , whose cumulants 
have a finite and unambiguous continuum limit HH mg. It is a naive discretization 
of the charge evolved with the Yang-Mills gradient flow. It is particularly appealing 
because its numerical evaluation is significantly cheaper than the one for the definition 
suggested by Neuberger’s fermions. Here we show that in the Yang-Mills theory the 
cumulants defined this way coincide, in the continuum limit, with those of the universal 
definition appearing in the anomalous chiral Ward Identities of QCD. By implementing 
the gradient-flow definition, we compute the topological susceptibility in the continuum 
limit with a precision 5 times better than the reference computation with the Neuberger’s 
definition M- We then determine the ratio of the forth cumulant over the second one 
in the continuum limit by keeping for the first time all systematics, especially finite 
volume effects, negligible with respect to the statistical errors. As a byproduct we also 
perform an interesting universality test at the permille level by comparing the values of 
the topological susceptibility at different flow times. 

2 Preliminaries in the continuum 

Starting from the ordinary fundamental gauge field 

-®/*| t=0 = > ( 2 * 1 ) 

where A p = A a JT a (see Appendix |A] for the generator conventions), the Yang-Mills 
gradient flow evolves the gauge field as a function of the flow time t > 0 by solving the 
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differential equation |12j 


d t B p 

- Bi/Gjyp + (M.qD , 


(2.2) 

Gpu 

= - d v B^ - i[BB v \ , 

D p = d p -i [B p , ■] , 

(2.3) 


with «o being the parameter which determines the gauge. Here we focus on the gradient- 
flow evolution of the topological charge density defined a^] 


q 


t 


1 

327 r 2 w tr 


G/ivGpa 


and of the corresponding topological charge 



(2.4) 


(2.5) 


where e pupa is the four-index totally antisymmetric tensor and the trace is over the color 
index. Under a generic variation 5B p of a given gauge held configuration, the topological 
charge density changes as 


V = dpwtp , 

^p — g^_2 \G pv 8B(j\ , 

(2.6) 

see for instance Ref. [15]. If we now 

specify 



5Bp = d t Bp St , 

(2.7) 

it is straightforward to show that 



d t q l = dpw'p , 

wp — g^_2 fr [G p jiy D a G acr \ , 

(2.8) 


where w p is a local dimension-5 gauge-invariant pseudovector held. This in turn implies 
that for a given gauge held configuration 


d t Q l = 0 , 


(2.9) 


an equation which reflects the topological nature of Q l . 

When q t (x) is inserted in a correlation function, Eq. (2.8) implies 

(q\x) Q(y )) = ( q t=0 (x ) Q(y)) + d p [ dt' {w* (x) 0(y)) (x / y) , 


( 2 . 10 ) 


where 0(y) is any finite (multi)local operator inserted at a physical distance from x. 
The l.h.s. of Eq. (2.10) is finite thanks to the fact that a gauge-invariant local composite 
held constructed with the gauge held evolved at positive how time is finite naE]. Since 


2 If not explicitly indicated, the superscript t on the quantities evolved with the gradient flow is 
always > 0. 
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there are no local composite fields of dimension d < 5 with the symmetry properties of 


Wp(x), the integrand on the r.h.s of Eq. (2.10) diverges at most logarithmically when 
t' —> 0. This implies that the quantity 


(<7 


,t= o, 


» °{y)) = Jim <?*(*) °(y)) (® + y ), 

t —>0 


( 2 . 11 ) 


is finite, i.e. the limit on the r.h.s exists for any finite operator 0(y). The Eq. (2.11) can 
be taken as the definition of q t=0 (x), i.e. the renormalized topological charge density 
operator at t = 0. The latter satisfies the proper singlet chiral Ward identities when 
fermions are included, see next section for an explicit derivation. It is worth noting that 


Eq. (2.10) implies that the small-t expansion of q t (x) is of the form 
{q\x) 0{y)) = {q t=0 (x) 0(y)) + 0(t) (x ^ y) . 


( 2 . 12 ) 


In the following we will be interested in the cumulants of the topological charge 


C* n = d 3 4 x i ... d A X 2 n-i{q\xi) ... q\x 2n )) c , 


(2.13) 


which, thanks to Eq. (2.9), are expected to be independent of the flow-time for t > 0 
with the limit t —> 0 which requires some care due to the possible appearance of short- 
distance singularities. It is the aim of the next section to address this question by using 
the lattice, a regularization where the theory can be non-perturbatively defined. We 
will supplement the theory with extra degenerate valence quarks of mass m, and we will 
consider the (integrated) correlator of a topological charge density with a chain made of 
scalar and pseudoscalar densities |7j defined as 


(q t =^0)P 5 i(z 1 )S 12 (z 2 )S 23 (z 3 )S u (z 4 )S^(z 5 )) 


(2.14) 


where S t j and Pij are the scalar and the pseudoscalar renormalized densities with flavor 
indices i and j. Power counting and the operator product expansion predict that there 
are no non-integrable short-distance singularities when the coordinates of two or more 
densities in (2.14) tend to coincide among themselves or with 0. When only one of 
the densities is close to q i=0 (0), the operator product expansion predicts the leading 
singularity to be 

q t=0 (x) Sij( 0) c(x) Pij{ 0) + 
where c(x) is a function which diverges as |x| -4 when |x| 


(2.15) 

0, and the dots indicate 
sub-leading contributions. An analogous expression is valid for the pseudoscalar density. 
Being the leading short-distance singularity in the product of fields q t=0 (x) Sij(0), its 


Wilson coefficient c(x ) can be computed in perturbation theory. By using Eq. (2.10), to 
all orders in perturbation theorjj^jwe can write 

(q t=0 (x)S ij (0)O(y)) = (q t (x)S ij (0) 0(y)) - d p [ dt , {w t '(x)S ij (0)O(y)) , (2.16) 


3 Since the function |as| 4 ln(a: 2 ) P is integrable for p > 1, the singularity needs to be determined 

only up to some finite order. 
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where again 0(y ) is any finite (multi)local operator inserted at a physical distance from 0 
and x. When t > 0, the first member on the r.h.s of Eq. (2.16) has no singularities when 
|x| —>• 0. If present, the singularity has to come from the second term, and therefore c{x) 
must be of the form 

c(x) = dpU p (x) (2.17) 


which does not contribute to the integral (over all coordinates) of the correlation function 

(2m 


3 Cumulants of the topological charge on the lattice 

On the lattice the Yang-Mills gradient-flow equation can be written as a first-order 
differential equation m 

d t Vp(x) = -gl{d x ,pS{V)}Vp{x), Vfi(x)\ t= o = Up(x) , (3.1) 

where the Wilson action S and the link differential operators d XtfX are defined in Ap¬ 
pendix [A] together with other conventions. The gauge field evolved at positive flow-time 
Vp(x) is smooth on the scale of the cut-off. When inserted at a physical distance, the 
gauge-invariant local composite fields constructed with the evolved gauge field are finite 
as they stand. Remarkably their universality class is determined only by their asymp¬ 
totic behavior in the classical continuum limit |121 [13] . At t > 0 any decent definition of 
the topological charge density is therefore finite. The same line of argumentation applies 
to the cumulants of the topological charge. At t > 0 short-distance singularities cannot 
arise because of the exponential damping of the high-frequency components of the fields 
enforced by the flow evolution. 

It remains to be shown, however, that the cumulants of the topological charge dis¬ 
tribution defined at t > 0 satisfy the proper singlet chiral Ward identities when fermions 
are included in the theory. To show this it is sufficient to work with a particular dis¬ 
cretization of the topological charge, and then appeal to the above mentioned universality 
argument for the other definitions. The GW discretizations have a privileged role since 
at t = 0 the lattice bare cumulants are finite, and they satisfy the singlet chiral Ward 
identities when fermions are included in the theory. 


3.1 Ginsparg Wilson definition of the charge density 

The definition of the topological charge density suggested by GW fermions is CS101EO 




a 

2 


tr 


75 D(x,x) 


(3.2) 


where we indicate it with a subscript N since, for concreteness, we take D(x,y) to be 
the Neuberger Dirac operator given in Appendix [A] in which each link variable Up(x) is 
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replaced by the corresponding evolved one V fl (x) when t > 0. Since there are no other 
operators of dimension d < 4 which are pseudoscalar and gauge-invariant, it holds that 


lim Z q (ql (0) °(x)) = finite 


(3.3) 


where Z q is a renormalization constant which is at most logarithmically divergent, while 
q^( 0) is finite as it stands. This in turn implies that 


lim Z a a 

a —>0 


E<? n(°) 


ql °( x )) = finite 


(3.4) 


since there are no short-distance singularities that contribute to the integrated corre¬ 
lation function because 5 * (0) is evolved at positive flow-time. By supplementing the 


theory with extra degenerate valence quarks of mass m, and by replacing in Eq. (3.4) 
the topological charge at t = 0 with its density-chain expression [7] we obtain 


a 




lit 


20 E' 

Zi,...,Z5 


'(<?n(0) -P 51 (~i) < 512 (^ 2 ) < 523 (^ 3 ) < 534 (^ 4 ) < 545 ( 2 : 5 )) , (3.5) 


where Sy and Py are the scalar and the pseudoscalar densities with flavor indices i and 
j. Written as in Eq. (3.5), power counting and the operator product expansion predict 


that there are no non-integrable short-distance singularities when the coordinates of two 


or more densities tend to coincide. The r.h.s of Eq. (3.5) is finite as it stands, and it 
converges to the continuum limit with a rate proportional to a 2 . This in turn implies 
that the limits on the l.h.s of Eqs. (3.3) and (3.4) are reached with the same rate if 


Z q is set to any fixed ( 50 -independent) value. Since in the classical continuum limit 
Neuberger’s definition in Eq. ( |3.2| ) has the same asymptotic behavior of the definition 
in Eq. ( | 2 ~T| ) PH EE], we may set Z q = 1 in which case 


lim {q N u (x) 0 L (y)) = (q u (x) 0(y)) (x / y) , 


a —>-0 


(3.6) 


where 0 L (y) is a discretization of the generic finite continuum operator 0{y). Once 
inserted in correlation functions at a physical distance from other (renormalized) fields, 
g(, = 0 (x) does not require any renormalization in the Yang-Mills theory. It is finite as 
it stands, and it satisfies the singlet Ward identities when fermions are included in the 
theory. It is interesting to note that Eqs. ( 2. 12| ) and (3.6) implies 

(q\x) 0{y)) = (q^°(x) 0 L (y)) + 0{a 2 ) + 0(t ) , (3.7) 

where in general discretization effects depend on t. We could have arrived to Eq. m b y 


following a procedure analogous to the one in the continuum, see Eqs. (2.8) -(2.12). To all 


orders in perturbation theory, or in general when Neuberger’s operator is differentiable 
with respect to the gauge field m, the change of the topological charge density with 


respect to the flow-time can be written, analogously to Eq. (2.8), as m f 21 ] (see also 


d t q N (x) = d*wl p (x) , (3.8) 

where w^ p (x) is a discretization of the dimension-5 gauge-invariant pseudovector oper¬ 
ator w 1 ' (x), and d* is the backward finite-difference operator. 
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3.2 Ginsparg Wilson definition of the charge cumulants 

The Neuberger’s definition of the topological charge is given by 

Ql = « 4 X/ ^0*0 , (3-9) 

X 

and its cumulants are defined as 


<n=« Sn - 4 £ {,'(li)...,'(l 2 „_ 1 )(,U0))c. (3.10) 

Xl,...,X2n-l 


For t = 0 the cumulants have an unambiguous universal continuum limit as they stand 
and, when fermions are included, they satisfy the proper singlet chiral Ward identities [5], 
017]. They are the proper quantities to be inserted in the Witten-Veneziano relations 
for the mass and scattering amplitudes of the rf meson in QCD POSED]- It is far from 
being obvious that coincide with those defined at positive flow-time, since the two 
definitions may differ by additional finite contributions from short-distance singularities. 

For the clarity of the presentation we start by focusing on the lowest cumulant, 
the topological susceptibility (7* 1 . At t = 0, by replacing one of the two g* =0 with its 
density-chain expression jT], we obtain 


°(9)^n °( x )) = - m5 a 20 'Y^(q t N °{0)P 51 (Z 1 )S12{Z2)S23(Z 3 )S 34: (Z4)S45(Z5)) ■ (3.11) 


Z1,...,Z5 


When the susceptibility is written in this form, the discussion toward the end of section [2] 


and in particular Eq. (2.17) guarantee that there are no contributions from short-distance 


singularities. This result, together with the fact that Z„ = 1, implies that 


lirn lima 4 ^(g*(x)4 °(0)) 

t -¥0 a —>0 z —' 

x 


a —>-0 z —' 

x 


(3.12) 


By replacing on the l.h.s g* °(0) with the evolved one, no further short-distance singu¬ 
larities are introduced and we arrive to the final result 


lim lima 4 ^(g*(x)g‘(0)) = lima 4 ^(g* °(x)g* °(0)) . 

0 a —>-0 z —' z ^ 


CL —^0 


(3.13) 


By replacing 2n — 1 of the charges in the n th cumulant with their density-chain 


definitions, the very same line of argumentation can be applied. The Eq. (3.13), together 
with the independence up to harmless discretization effects of Cj. n from the flow-time 
for t > 0 |12j . implies that the continuum limit of n coincides with the one of <7^^. 
The cumulants of the topological charge distribution defined at t > 0 thus satisfy the 
proper singlet chiral Ward identities when fermions are included ppm- They are the 
proper quantities to be inserted in the Witten-Veneziano relations for the mass and 
scattering amplitudes of the r/ meson in QCD PIPED]. 
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3.3 Universality at positive flow-time 

For t > 0 different lattice definitions of the topological charge density belong to the same 
universality class if they share the same asymptotic behavior in the classical continuum 
limit Hama. In the rest of this paper we are interested in the naive definition of the 
topological charge density defined a^] 

Q\ x ) = ^2 G%( x )Gpa( x ) , (3-14) 

where the field strength tensor G^flx) is defined as [23] 

G^(x) = tr l(Q^(x) - Q„u{x)) T a ] , (3.15) 

with 


Quips') 


V^{x)V u (x + ajl)V^(x + av)V$(x) + 

V u (x)V^(x - afi + av)Vj(x — affjV^x — afl ) + 

V^(x — afl)Vf (x — afi — az>)V^(x — afi — av)V v (x — av) + (3.16) 

VJ(x — au)V^ t (x — av)V y {x + afl — av)V^(x) . 


In the Yang-Mills theory q t °(x) requires a multiplicative renormalization constant]^ 
when inserted in correlation functions at a physical distance from other operators 


The cumulants of the corresponding topological charge Q t=u = uJJ x Q t=U (x), defined 
analogously to Eq. (3.10), have additional ultraviolet power-divergent singularities, and 
they do not have a well defined continuum limit. 

The density q t (x ) in Eq. (3.14) shares with q^(x) the same asymptotic behavior in 
the classical continuum limit mm- Since for t > 0 short-distance singularities cannot 
arise, C £ n and C\ tend to the same continuum limit. The results in the previous section 
then imply that the continuum limit of the naive definition of C*, at positive flow-time, 
coincides with the universal definition which satisfies the chiral Ward identities when 
fermions are added E El 0. It is interesting to note, however, that at fixed lattice 
spacing there can be quite some differences. For instance, the topological susceptibility 
defined at t > 0 with the naive definition is not guaranteed to go to zero in the chiral 
limit at finite lattice spacing in presence of fermions [25]. 


4 Numerical setup 


For the numerical computation we discretize the SU(3) Yang-Mills theory with the 
standard Wilson plaquette action on a finite four-dimensional lattice with spacing a, 


4 We use the same notation for the naive definition of the field strength tensor, of the topological 
charge and of its density on the lattice and in the continuum, since any ambiguity is resolved from the 
context. 

5 This renormalization constant can be fixed by enforcing the analogous of Eqs. (3.31 and (3.6 1 . 
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Lattice 

P 

to/a 2 

L/a 

L [fm] 

a [fm] 

A^conf 

nit 

eerr 

q2err 

q4err 

Ai 

5.96 

2.79 

10 

1.0 

0.102 

36 000 

30 

0.19 

0.0005 

0.0024 

Bi 



12 

1.2 


144000 


0.45 


0.005 

Ci 



13 

1.3 


280 000 


0.42 


0.0068 

Di 



14 

1.4 


505 000 


0.74 


0.01 

Ei 



15 

1.5 


880 000 


0.89 


0.012 

Fi 



16 

1.6 


1440 000 


1.04 


0.015 

b 2 

6.05 

3.78 

14 

1.2 

0.087 

144000 

60 

0.31 

0.0005 

0.005 

d 2 



17 

1.5 


144000 


0.045 


0.01 

b 3 

6.13 

4.87 

16 

1.2 

0.077 

144000 

90 

0.25 

0.0005 

0.005 

d 3 



19 

1.5 


144000 


0.058 


0.01 

Bi 

6.21 

6.20 

18 

1.2 

0.068 

144000 

250 

0.20 

0.0005 

0.005 

T>4 



21 

1.4 


144000 


0.042 


0.01 


Table 1 : Overview of the ensembles and statistics used in this study. For each lattice we 
give the label, /3 = 6 /« 7 q, the reference scale to/a 2 , the spatial extent of the lattice, the 
lattice spacing, the number !V con f of independent configurations generated, the number 
of sweeps nit required to space them, and the tolerances eerr, q2err and q4err on the 
primary observables considered (see main text). 


with the same L/a size in all four space-time directions, and with periodic boundary 
conditions imposed on the gauge fields, see Appendix [A] for details. The basic Monte 
Carlo update of each link variable implements the Cabibbo-Marinari scheme [26] . by 
sweeping the full lattice with one heatbath update followed by L/(2a) sweeps of over¬ 
relaxation updates. 


4.1 Ensembles generated 


We have simulated three series of lattices in order to estimate and remove the systematic 
effects due to the finiteness of the lattice spacing and volume, see Table [l] for details. 
In the first series {A \, B \,..., F\} the inverse coupling /3 = 6 /^q is kept fixed so that 
the lattice spacing is approximative^ 0.1 fm, while the physical volume increases from 
(1.0 fm ) 4 to (1.6 fm) 4 . The number N con f of independent gauge configurations generated 
scales with L 8 to ensure that the relative statistical error on R, the ratio of the fourth 


over the second cumulant of the topological charge distribution see Eq. (4.2), is always 
at the 10 % level m ■ In the second series {B i,..., B 4 } the physical volume is kept 
approximative^ fixed, while the spacing is decreased down to 0.068 fm. The volume is 
always (1.2 fm ) 4 to guarantee that finite-size effects 011 R are within the statistical errors, 
while the computational cost remains affordable. In the third series {Di,..., D 4 } is 
again the physical volume which is kept approximatively fixed, always at least (1.4 fm) 4 , 
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to guarantees that finite-size effects in the reference scale to and in the topological 


susceptibility y, see Eq. (4.2), are within their (smaller) statistical errors. In both cases 


the measurements at the four lattice spacings are used to estimate discretization effects 
in the observables, and to extrapolate them away in the continuum limit. 


4.2 Computation of the observables 

The primary observables that we have computed on each configuration at t > 0 are the 
energy density 




and the topological charge Q defined as in section |3.3 
interested in are 

m 2 ) 


{E l 


X* = 


$(*) > ( 4 - 1 ) 

The quantum averages we are 


V 


R t = 


([Q 


1 12 \ 


(4.2) 


To numerically integrate the Yang-Mills gradient flow we have implemented a fourth or¬ 
der Runge-Kutta-Munthe-Kaas (RKMK) method [281 l29l, 130 J. It is a structure-preserving 
Runge-Kutta (RK) integrator, designed to exactly preserve the Lie group structure of 
the gradient flow equation, see Appendix [B] for details. On each lattice the field has 
been evolved approximatively up to t = 1.2 to, where to is the reference flow-time value 


defined below. The observables in Eq. (4.2) have been computed with a flow-time reso¬ 


lution of 0.08a 2 or smaller. The numerical integration of the flow equation introduces a 
systematic error in the gauge field val¬ 
ues at positive flow-time t, and thus in 
each observable. In our case, at asymp¬ 
totically small values of the RK step 
size, it is proportional to e 4 . There are, 
however, large fluctuations in the pre- 
factor among the various gauge configu¬ 
rations, see Figure [lj A reliable estimate 
of this systematics is achieved by moni¬ 
toring the error configuration by config¬ 
uration, and occasionally adapt the step 
e. To do so we integrate the flow equa¬ 
tion two times with steps e and 2e, where 
in our case e = 0.08a 2 . Denoting with 
E and the basic observables E and 



Q respectively computed on the j th field 
configuration evolved with step size e, at 
small enough e the error is given by 


#cfg 

Figure 1: History plot of the systematic error 
5Q 2 4 2e ) at t ~ to f° r the first 2000 configura¬ 
tions of D 4 . The red line indicates the bound 
q2err = 0.0005 of the systematic error en¬ 
forced on all configurations. 
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5E (2e) = 


^>( e ) _^(2c) 


5Qf e) = 


Qf - Q? e) 


(4.3) 


with both observables obviously measured at the same flow-time. By applying linear 
propagation, the error on the average over all configurations is bounded by 

5Ef e) < max (<5^f e) ) , <5Q 2 ’ (2e) < max (|2Q„ I^Q^) , ^0 4,(2e) < max (|4 Q^SQ^ 


SR^ < i max ^max ® + Q 2 ^ ^ max (l • (4.4) 

At run-time, for each configuration and each flow-time the systematic errors of the 
observables E, Q 2 and Q 4 are compared with the given tolerances eerr, q2err and 
q4err respectively. If one of the tests fails, the flow evolution is re-computed for that 
configuration with a new step size e' = (l/ 2 )e and new observables data, along with old 
e = 2 e' data, are used to estimate the systematic errors and compare them with the 
tolerances. If the test fails again, the field is evolved with e" = (l/2)e / , and so on. This 
ensures that 


SE (2e) 


j — 


< eerr 


<5Q 2 ’ (2e) < q2 


err 


dQ 4 ’*' 2 ^ < q4err . 


(4.5) 


The parameters eerr, q2err and q4err are chosen as a function of the target statistical 
error on the corresponding observables. If we set the upper limit for the systematic 
error to be roughly 10 times smaller than the statistical one, this condition is readily 
translated into a limit for eerr, q2err and q4err, see Table [l] for the values chosen for 


each lattice. The Eqs. (4.5) put bounds 


on the systematic errors for the coarser 
evolution, but the data evolved with the 
finer step size e are actually those used in 
the final analysis. This choice is rather 
conservative in our case, being the actual 
error more than one order of magnitude 
smaller. For the quantity E t the actual 
error turns out to be more than two or¬ 
ders of magnitude smaller with respect 
to the bound in Eq. (4.5), see Figure [ 8 ] 
in Appendix [B] We have therefore chosen 
larger values for eerr with respect to one 


given by the bound in Eq. (4.5). 


4.3 Autocorrelation times 



(a 2 /to)-h 2 

Figure 2: The integrated autocorrelation 
times ri nt of the primary observables as a 
function of (a 2 /to) -1 ^ 2 - 


To measure the autocorrelation time of the various observables, we perform a dedicated 
run for each lattice {-Bi,..., B 4 } where the gauge field configurations are separated by 
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Lattice 

-^conf 

t/a 2 

T . E 
' lilt 

r Q 

int 

r Q2 

'int 

r Q4 

1 int 

B\ a 

36 x 1000 

3.36 

7.0(5) 

9.1(7) 

5.2(3) 

4.39(25) 

L>2a 

36 x 1000 

4.64 

7.9(6) 

17.4(18) 

7.9(6) 

6.3(4) 

L>3 a 

36 x 1000 

6.08 

12.4(11) 

30(4) 

12.4(11) 

9.2(7) 

B^a 

36 x 1000 

7.68 

17.2(18) 

73(13) 

31(4) 

22.1(25) 


Table 2: Integrated autocorrelation times of the various observables in units of a single 
sweep of the update algorithm. They have been measured on dedicated runs made of 36 
series of 1000 sweeps each. 


a single iteration of the update algorithm. Each series is replicated 36 times to increase 
statistical accuracy. The integrated autocorrelation times r; n t of the observables E, Q , 
Q 2 , and Q 4 , estimated as in Ref. |3T| . are reported in Table [ 2 ] In the range of /3 values 
considered, Q has the largest autocorrelation time which increases rapidly toward the 
continuum limit |32| . To ensure that the measurements in the main runs are statistically 
independent, we have spaced them by nit sweeps of the lattice, see Table [lj 

5 Physics results 

A first analysis of the data reveals the effectiveness of the gradient flow in splitting the 
field space of the lattice theory into different topological sectors. In Figure [3] we plot the 
histograms of the topological charge Q measured at different flow-times on the lattice 
D4. In the plot on the top-left corner, the topological charge distribution at t = 0 is a 
smooth function over non-integer values. By increasing the flow-time, the configurations 
with charge close to integers become more and more probable. The spikes in the bottom- 
right plot turn out to be slightly shifted towards zero with respect to the integer values 
due to discretization effects. On the other lattices similar histograms are obtained. 


5.1 Scale setting 

The reference flow-time to is defined through the implicit equation [IT2] 


t 2 (-E f )| ( _. — 0.3 . 


(5.1) 


In the region of interest t 2 (-E 4 ) grows approximatively as a linear function of t. Since we 
have computed (-E*) at flow-times spaced by finite steps, we have solved equation (5.1) 
by interpolating linearly the two data points closest to to- The results are reported 
in Table [3j with the systematic error due to the interpolation being negligible. By 
comparing the values of to/a 2 obtained on the lattices {A \,..., Ti}, finite-size effects 
are not visible at the level of 0.1 permille in the statistical precision for L > 1.4 fin. 
We thus fix the lattice spacing at all values of /3 from to/a 2 determined on the lattices 
{ D\ ,..., D 4 }. 
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Figure 3: Histograms of the topological charge distribution measured on the lattice D 4 
at different flow-times. 

In Table [ 3 ] the values of to/ r o> w h ere r o is the Sommer scale computed in [33], 
are also reported. As shown in Figure [4| discretization effects in this ratio are indeed 
negligible with respect to the statistical errors dominated by the 0.3-0.6% error on r^/a. 
By extrapolating the results linearly in a?/to, we obtain in the continuum limit 

V^to/ro = 0.941(7), (5.2) 

which corresponds to to /tq = 0.1108(17). To express to in physical units, we supplement 
the theory with quenched quarks. The value of F^r 0 = 0.293(7) from Ref. [34] together 
with Fk = 109.6 MeV leads tc(2] 

t 0 = (0.176(4) fm) 2 , (5.3) 

the error being dominated by the one on FkT'o. 

e Note that we use an updated determination for the physical value of Fk with respect to Ref. |3!j . 
The change on Fkt 0 induced by the new tuning of the strange quark mass is negligible with respect to 
the statistical error quoted. 
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Lattice 

t 0 /a 2 

*oAo 

Ai 

2.995(4) 

0.1195(9) 

Bi 

2.7984(9) 

0.1117(9) 

Ci 

2.7908(5) 

0.1114(9) 

D\ 

2.7889(3) 

0.1113(9) 

Ei 

2.788 92(23) 

0.1113(9) 

F\ 

2.788 67(16) 

0.1113(9) 


Lattice 

to/a 2 

to/rl 

b 2 

3.7960(12) 

0.1114(9) 

Bo 

4.8855(15) 

0.1113(10) 

Bi 

6.2191(20) 

0.1115(11) 

d 2 

3.7825(8) 

0.1110(9) 

Do 

4.8722(11) 

0.1110(10) 

d 4 

6.1957(14) 

0.1111(11) 


Table 3: Results for the reference flow-time to/® 1 2 and the ratio to/r q. The error on the 
latter is dominated by the 0.3-0.6% relative error on ro/a quoted in |33] , 


5.2 Topological susceptibility 


The full set of results for the topological charge moments and cumulants are given in 
Table [ 4 ] They are computccj/] at the reference flow-time to by linearly interpolating the 
numerical data as described in the previous section. 

In Figure [5] we show the values of the 
topological susceptibility y = (Q 2 )/V 
from the lattices {A \,..., F\\ as a func¬ 
tion of the linear extension of the lat¬ 
tice. For L > 1.4 fm, finite-size effects 
turn out to be below our target statisti¬ 
cal error of approximatively 0.5%. The 
continuum value of t/y can thus be ob¬ 
tained by extrapolating the results from 
the lattices {D 1 ,..., D 4 }, see left plot of 
Figure [6] The Symanzik effective theory 
analysis predicts discretization errors to 
start at 0(o 2 ), and indeed the four data 
points are compatible with a linear be¬ 
havior in a 2 . A linear fit of all of them 
gives as intercept t^y = 6.75(4) x 10~ 4 
with a significance of y 2 /dof = 1.26. A 
quadratic fit gives t(jy = 6.49(18) • 10 -4 
with y 2 /dof = 0.38, and with a coeffi¬ 
cient of the quadratic term compatible 

with zero within the statistical errors. By restricting the linear fit to the three points at 



a?/to 

Figure 4: Continuum limit extrapola¬ 

tion of y/8to /ro computed 011 the lattices 
{D 1 ,..., D 4 }. The errors are dominated by 
the 0.3-0.6% relative error on ro/a quoted 
in 1331. 


1 Unless explicitly indicated, the gradient flow-time at which the topological quantities are computed 

throughout this and the next sections is t = to- 
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Lattice 

(Q 2 ) 

(Q 4 ) 

{Q A )c 

R 

Ai 

0.701(6) 

1.75(4) 

0.273(20) 

0.39(3) 

B i 

1.617(6) 

8.15(7) 

0.30(4) 

0.187(24) 

Ci 

2.244(6) 

15.50(10) 

0.40(5) 

0.177(23) 

Di 

3.028(6) 

28.14(14) 

0.63(7) 

0.209(23) 

Ei 

3.982(6) 

48.38(18) 

0.81(9) 

0.202(23) 

Fi 

5.167(6) 

80.90(22) 

0.81(11) 

0.157(22) 

b 2 

1.699(7) 

9.07(9) 

0.41(5) 

0.24(3) 

d 2 

3.686(14) 

41.6(4) 

0.83(19) 

0.22(5) 

b 3 

1.750(7) 

9.58(9) 

0.39(5) 

0.22(3) 

d 3 

3.523(13) 

37.8(3) 

0.56(17) 

0.16(5) 

b a 

1.741(7) 

9.44(9) 

0.35(5) 

0.20(3) 

d 4 

3.266(12) 

32.7(3) 

0.68(15) 

0.21(5) 


Table 4: Results for the various topological observables measured at flow time to on all 
lattices simulated. 

the finer lattice spacings, we obtain 

to X = 6-67(7) x lO" 4 , (5.4) 

with y 2 /dof = 0.88, which is our best result for this quantity. It is five times more 
precise than the determination which uses the Neuberger’s definition of the topological 
charge m- 

The cumulants of the topological 
charge are expected to be f-independent 
in the continuum limit. In the right 
plot of Figure [6] we show the topological 
susceptibility computed at various flow- 
times normalized to its value at to- The 
data points have statistical errors which 
range from 0.1 to 1 permille due to the 
correlation between the numerator and 
the denominator. At finite lattice spac¬ 
ing discretization effects are clearly vis¬ 
ible, and they depend on t. When each 
set of data is extrapolated to the con¬ 
tinuum limit with a quadratic function 
in a? /to, the intercepts are all compat¬ 
ible with 1 within the statistical errors 
which, depending on t, range from 0.5 to 



Figure 5: Values of a 4 y as a function of L for 
the series {Ai, ... ,Fi}. 
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a?/to a 2 /to 


Figure 6: Right: the dimensionless quantity x as a function of a 2 /to, and its extrapo¬ 
lation to the continuum limit. Left: the ratio xVx (errors are smaller than symbols) as 
a function of a 2 /to for several values of t, and its extrapolation to the continuum limit. 


5 permille. We can also compare our result in Eq. (5.4) with the one obtained almost 
10 years ago with the Neuberger’s definition of the topological charge M- If we use 
Eqs. (5.2) and (5.4), we obtain 


4x = 0.0544(18), 


(5.5) 


which differs by less than 1.5 standard deviationsj^jfrom the result in Eq. (11) of Ref. [[T4] , 
It is interesting to note that after ten years from the first computation of x i n the 
continuum limit m, we moved from an unsolved problem to a universality test at the 
permille levej^J 


By using the result in Eq. (5.3), the value of x hi physical units is given bj 


10 


X = (180.5(5)(43) MeV) 4 


(5.6) 


where the first error is statistical from Eq. (5.4), while the second is the one from the 
uncertainty in the scale in Eq. (5.3). If we use the physical value of to determined in 
QCD with Nf = 2 and Nf = 2 + 1 flavours (3?| l38j . we obtain a value of x hi physical 
units which differs (downwards) by 10-20% per linear dimension. This is the size of 
the ambiguity which is expected when results of the Yang-Mills theory are expressed in 
physical units. 


8 This value takes into account the fact that the same determination of ro is used in the two compu¬ 
tations. 

9 A first test of universality for x was already presented in Ref. [35] with statistical errors more than 
one order of magnitude larger than those obtained here. Results with similar large statistical errors 
were recently obtained in Ref. (36) . 

10 Note that in Ref. El. Fk = 113.1 MeV was used to set the scale in physical units. 
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Figure 7: Left: values of R at flow-time to versus L for the series {A \...., Fi}. Right: 
the quantity R as a function of a 2 /to and its extrapolation to the continuum limit; the 
dotted blue hue is the dilute instanton gas model prediction R = 1. 


5.3 The ratio R 

The values of R = (Q 4 ) c /( Q 2 ) from the lattices {A\, ..., Fi} are shown in the left plot of 
Figure [7] as a function of L. Since our target statistical error is approximatively 10%, a 
linear extension of L > 1.2 fm is enough for finite-size effects to be within errors. Given 
the increase with L 8 of the computational cost of R, we have chosen to determine its 
continuum limit by extrapolating the data from the lattices {B i,..., B 4 }, see left plot of 
Figure [7] Also in this case the Symanzik effective theory analysis predicts discretization 
errors to start at 0(a 2 ), and indeed the four data points are compatible with a linear 
behavior in a 2 . A fit to a constant of all of them gives R = 0.210(13) with a significance 
of y 2 /dof = 0.83. A linear fit in a 2 /to gives 


R = 0.233(45), 


(5.7) 


which is our best result for this quantity. The significance of the fit is \ 2 /dof = 1.1, and 
the slope is compatible with zero. 

The value in Eq. (5.7) is compatible with the one obtained with the Neuberger’s 
definition in Ref. m , albeit with an error 2.5 times smaller. It is also relevant to 
note that a systematic study of finite-size effects was not carried out in Ref. |27j . and 
finite-size effects were estimated and added to the final error. 


6 Conclusions 

The 0-dependence of the vacuum energy, or equivalently the functional form of the topo¬ 
logical charge distribution, is a distinctive feature of the ensemble of gauge configurations 
that dominate the path integral of a Yang-Mills theory. The value of R = 0.233(45) in 
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Eq. (5.7) rules out the 0-behavior predicted by the dilute instanton gas model. Its large 
distance from 1 implies that, in the ensemble of gauge configurations that dominate the 
path integral, the fluctuations of the topological charge are of quantum non-perturbative 
nature. The large N c expansion does not provide a sharp prediction for R. Its small 
value, however, is compatible with being a quantity suppressed as 1/N% in the limit of 
large number of colors N c . The value of R found here is related via the Witten- Veneziano 
mechanism to the leading anomalous contribution to the rf—rf elastic scattering ampli¬ 
tude in QCD. It is one of the low-energy constants which enter the effective theory of 
QCD when its Green functions are expanded simultaneously in powers of momenta, 
quark masses and 1 /N c . 

The Yang-Mills gradient flow is an extremely powerful tool for studying the topo¬ 
logical properties of the theory. It provides a reference scale and a sensible definition 
of the topological charge which are cheap to be computed numerically. With a modest 
numerical effort by today standards, it allowed us to compute the dimensionless ratio 
t^x = 6.67(7) x 10 -4 with a relative error of roughly 1% in the continuum limit, i.e. 
five times smaller than the one of the previous reference computation with the Neu- 
berger’s definition. The Yang-Mills gradient flow is clearly an interesting tool to study 
the topological properties of the Yang-Mills vacuum as a function of N c . 

As proven in this paper, in the continuum limit the cumulants of topological charge 
defined by the Yang-Mills gradient flow coincide with those of the universal definition 
appearing in the chiral Ward identities. This in turn implies that this definition of 
the topological charge is the correct one for studying the 0-dependence of the vacuum 
of QCD at zero and non-zero temperature. If computed in thermal (full) QCD, its 
cumulants can be directly related, for instance, to the axion dynamics without further 
renormalization. 
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A Definition and conventions 

The Lie algebra of the SU(3) group may be identified with the linear space of all hermi- 
tian traceless 3x3 matrices. In the basis T a , a = 1... 8, with 

tr[T a ] = 0 , T at = T a , (A.l) 
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the elements of the algebra are linear combinations of them with real coefficients. The 
structure constants f abc in the commutator relation 

[T a , T b ] = i f abc T c (A. 2) 


are real and totally anti-symmetric in the indices if the normalization condition 

1 


tr 


rj-iCLrpb 


= -5 ab 
2 


(A.3) 


is imposed. 

For the SU(3) Yang-Mills theory the standard Wilson plaquette action is given by 


S i U ] = 1 - ^Retrjt/^x)} 


X LL,1S 


(A.4) 


where the trace is over the color index, j3 = Q/g$ with go the bare coupling constant, a 
is the lattice spacing, and the plaquette is defined as a function of the gauge links U^(x) 
as 

Ufax) = U^x) U v {x + ajl) Ufa + av) Ufa) , (A.5) 

with n, v = 0,..., 3, fi is the unit vector along the direction g and x is the space-time 
coordinate. 

The Neuberger-Dirac operator is defined as |Tj 

D = 3 {1 + 75 sign (H)}, 

a 

H = 75 ( ciD w - 1 - s) , a = fa— , (A. 6 ) 

1 + s 

where s is a real parameter in the range |s| < 1, and and D w is the the Wilson-Dirac 
operator. It is defined as 

Av = \ { 7 m(V; + V„) - aV* V M } , (A.7) 

where 

V M /(x) = - {Ufa)f(x + ap) - f(x)} , (A. 8 ) 

a 

V*/(.t) = ^fa(x)-Ufa - afi)f(x - afa (A.9) 

are the gauge-covariant forward and backward difference operators. The Neuberger- 
Dirac operator satisfies the GW relation 


7 5D + -D75 = 0D75D . 


(A.10) 
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The link differential operators acting on functions f(U) of the gauge field are 


d 


d a x J{U) = f(e~ isX U) 


X,fl 


d s 


x ( \ = ) Ta if x >**) 

s= o ^ I 0 otherwise 


While these depend on the choice of the generators T a , the combination 

d x ,,f(U) = T-d^fiU) 

can be shown to be basis-independent. 

B Runge—Kutta—Munthe-Kaas integrators 

Consider an ordinary differential equation 

y = f{y)y, y(fi) = yo, 


(A. 11 ) 


(A.12) 


(B.l) 


where y € G for some Lie group G and f{y ): G —> 0 , with g being the Lie algebra of G. 
Runge-Kutta-Munthe-Kaas methods |28l 1291 30] are structiLre-preserving Runge-Kutta 
methods designed to integrate numerically these equations on the group manifold, for a 
general introduction see Ref. [39]. The starting point is to write the solution of (B.l) as 


y(t) = exp {f (t)} y( 0 ) , 

and then solve the ordinary differential equation 

v = dexp " 1 {/(?/)} , r( 0 ) = 0 , 

where dexp ” 1 has the series expansion 


°° ~D 11 

exp" 1 = J[ ad v = 1 + 2 '] + 'll + 


(B.2) 


(B.3) 


(B.4) 


k =0 


with being the Bernoulli numbers, and ad^ = [w, •] the adjoint action. Since v(t) takes 


values in the Lie algebra, the differential equation (B.3) can be numerically integrated 


using an ordinary RK method. No extra conditions are needed, and any RK method of 
a given order can be used as a base for a RKMK method of the same order. The only 
complication is given by the operator dexp. 


v 1 , which can be substituted with its series 


expansion in Eq. (B.4) suitably truncated according to the order of the method. The 
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RKMK method of q th order with s stages is given by 


for i = 1, 2,..., s : 

S 

Ui — ^ ^ / ■] 

3 =1 

h = hf{exp(ui)y 0 } 

ki = dexpinv(rii, ki , q ) 
s 

v = ^ 2 b jk 

3 =1 

2/1 = exp{r}y 0 


where dexpinv(rt, u, g) is the truncated series 


9-1 

dexpinv(u, v, q) = 

k =0 



k 

U 


(B.5) 


(B.6) 


and cijj, are the coefhcients of g th order s-stages RK method. The fourth order 
RKMK method that we implemented is obtained starting from the very common 4 th 
order 4-stages RK method with coefficients, arranged in a Butcher tableau, 


0 


1 

1 

2 

2 

1 

2 

0 \ 

1 

0 0 1 


1111 


6 3 3 6 


(B.7) 


introduced by Kutta himself. At a first sight this method entails the computation of 
six different commutators of ki structures. However, it is possible to reduce the number 
of independent commutators needed to only two. As explained in Ref. m, this is 
due to the fact that, whereas the ki are in general O(h), some combinations of them are 
higher order in h and so the corresponding commutators can be neglected. The resulting 
integration algorithm is 


u\ = 0 , ki = hf{exp(ui)y 0 } 

U2 = ~k \, 

u 3 = ^k 2 + ^[ki,k 2 ] , 
u 4 = fa , 

v = \ki + \k 2 + \fa + \fa - — [fa, fa] , 

o 6 6 o 12 

2/1 = exp(u)y 0 . 


20 




Alternative RK methods for integrating 


(B.l) are given by the Crouch-Grossman 
integrators HD 02]. They are a spe¬ 
cial case of so-called commutator-free Lie 
group methods |43| . The third order al¬ 
gorithm described in Ref. [T2] belongs 
to this class. The conditions which the 
coefficients need to satisfy, order by or¬ 
der, are computable up to arbitrary or¬ 
der @D- They are given by the order con¬ 
ditions for a classical RK method, plus 
specific extra conditions. At fourth or¬ 
der, however, we did not find a coefficient 
scheme with the useful properties of the 
Lrischer’s integrator in terms of exponen¬ 
tial reusing. 


B.l Application to 
Mills gradient flow 


the Yang 


nr - 
io 2 - 
io° - 
K| io -2 - 
io -4 - 
io -6 - 
10 -8 — 


t 


t euler 
+ rk 2 
1 rk 3 
1 rk 4 mk 


nr 


10 ” 


Figure 8: Comparison of the numerical in¬ 
tegration methods. The systematic error 
6 E t ^ e \ where t = 3.2a 2 , is estimated from 
100 configurations of the lattice B\ by tak¬ 
ing the difference between E t evolved with 
step size e, and E t ( £ / 2 \ evolved with e/2. 


The Yang-Mills gradient flow equation (3.1) can be written as an ordinary first-order 
autonomous differential equation 


V(t) = Z[V(t)]V(t), V(0) = v 0 , 


(B.9) 


where 

Z[V{t)} = -g 2 0 {d x ,„s[v(t)}} , (B. 10 ) 

and the link differential operators are defined in Eq. ( |A.11| ). The fourth order RKMK 
method in (B.8) reads 


Wi = V(t) , Zi = eZ[Wi\ , 

W 2 = exp |^i | V(t) , 

W 3 = e W S^Z 2 + ^[Z 1 ,Z 2 ]\v(t) , (B.ll) 

W 4 = exp {Z3} V (t) , 

V(t + a 2 e) = exp |-Zi + -Z 2 + -Z 3 + -Z 4 — — [.Zf, Z 4 ] \ V(t) . 


This method computes four times the force held Z[Wi ] and four times the Lie group 
exponential. The commutators are economically implemented exploiting structure con¬ 
stants of 0 . Each iteration needs space in memory for one auxiliary gauge held and 
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three Z % fields. Gauge fields are stored in memory with a full 3x3 complex matrix, 
which has 18 real components, for each link. A Zj field is an element of su(3), which is 
a 8 dimensional linear space, for each link. Thus, the method (B.ll) requires space for 
(18 + 3 x 8 ) x 4V floating point numbers. Each exponential of a 0 -valued combination of 
Z[Wi\ reduce to 4V su(3) matrices exponentials, which can be computed economically 
exploiting the Cayley-Hamilton theorem as described in [ 35] . 

In the left plot of Figure [ 8 ] the RKMK method is compared to lower-order Runge- 
Kutta methods, such as the third order method rk3 found in [12]. The comparison is 
done averaging over 100 configurations at /3 = 5.96 on a 12 4 lattice evolved at t = 3.2 a 2 . 
The rk4mk algorithm scales correctly as a fourth-order method. However, the pre-factor 
appears to be larger, thus the new method is more precise with respect to rk3 in Ref. m 
for e < 0 . 1 . 
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